Orbital Magnetism in Three-Dimensional Quantum Dots 
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We study orbital magnetism in a three-dimensional (3D) quantum dot with a parabolic confining 
potential. We calculate the free energy of the system as a function of the magnetic field and the 
temperature. By this, we show that the temperature-field plane can be classified into three regions 
in terms of the characteristic behavior of the magnetization: the Landau diamagnetism, de Haas- van 
Alphen oscillation and mesoscopic fiuctuation of magnetization. We also calculate numerically the 
magnetization of the system and then the current density distribution. As for the oscillation of 
the magnetization when the field is varied, the 3D quantum dot shows a longer period than a 2D 
quantum dot which contains the same number of electrons. A large paramagnetism appears at low 
temperatures when the magnetic field is very weak. 
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PACS numbers: 



Orbital magnetism of a bulk system appears as the 
Landau diamagnetism (LD). In the past. Landau's 
work has provoked discussion regarding his treatment 
of boundary |^ and the effects of boundary on the 
orbital magnetism have been studied by many peo- 
ple. [|[ ^ Recently, two-dimensional (2D) quantum dots 
have become experimentally available with the capabil- 
ity of controlling the dots' size, shape and number of 
electrons. ^ Stimulated by experimental develop- 
ment, many theorists have been studying on 2D quan- 
tum dots with parabolic confining potentials. In partic- 
ular, Ishikawa and Fukuyama (IF) |^ have claimed that 
the temperature-field (T-B) plane can be classified into 
three regions in terms of the characteristic behavior of the 
magnetization M{T,B): the LD, de Haas- van Alphen 
(dHvA) oscillation and mesoscopic fluctuation (MP) of 
magnetization. Such magnetic behavior is mainly gov- 
erned by the energy levels near the Fermi level of the 
system. In the case of 3D quantum dots, the additional 
degree of freedom makes the density of states larger than 
that of 2D quantum dots. Hence, it is expected that the 
temperature affects the magnetization more easily, which 
changes the magnetic behavior to some extent and pos- 
sibly invalidates the classification. 

In this letter, we present our study on the orbital mag- 
netism of 3D quantum dots. We consider noninteracting 
electrons confined in a parabolic potential and calculate 
analytically the free energy of the system. Thereby, 
we show that the temperature-field plane can still be 
classified into the above three regions in terms of the 
characteristic behavior of the magnetization, M{T, B) — 
—d^l/dB. In order to confirm this classification, we cal- 
culate numerically the magnetization of the system and 
then the current density distribution. In addition to the 
consistency between the analytical and numerical calcu- 
lations, we report some remarkable findings of our nu- 
merical calculation. Hereafter, we set the fundamental 
constants as h = c = = ^■ 

We consider spinless electrons of effective mass m and 
charge — e confined in a spherical quantum dot with a 



parabolic confining potential in a uniform magnetic field 
B. The total Hamiltonian H of the system is given by 
the sum of single electron Hamiltonians as 



H = 



N 

E{2^HV. + eA(r,)f +fc.|r2}, (1) 



where N is the total number of electrons, and A is a 
vector potential satisfying B = rotA. The last term 
represents the confining potential, where W3 parametrizes 
the strength of the isotropic confinement. 

The eigenstates of the single electron Hamiltonian 
comprise the Fock-Darwin states in the x-y plane and 
the eigenstates of a harmonic oscillator in the z-direction 
parallel to B. In the cylindrical coordinate (p, 0, z), the 
eigenfunction ip is given by 

i^naA^) = (27r)~5 exp(ia0) Rna{p) Z,,{z), (2) 

where n, a, v are integers and n, v > 0. (n, a, ly) consists 
of a set of quantum numbers. Rna{p) is the radial wave 
function of the Fock-Darwin state specified by {n, a) and 
Z^{z) is the i^-th eigenfunction of the harmonic oscillator. 
The eigenenergy Enaf is given by 



Enav = {n+^)uJi + \a\ UJ2 + (i^ +5) ^3, 



W2(±) = t^± = 5(^1 ± U>b), 



(3) 

(4) 
(5) 



where cub denotes the cyclotron frequency, lub ~ eB/m. 
In eq. (^), the characteristic frequency uj2 — ^2(0") de- 
pends on the sign a {— ±) of the angular momentum a. 
Hereafter, we set uj^ as the energy unit, though we keep 
the symbol lj^ in the several equations below for clarity. 

The magnetization AI of the entire system is given by 
M ~ —{dVt/dB)Tf_i, Q where T is the temperature, /j, 
is the chemical potential, and Vl is the free energy of the 
system, n = -TE„a.ln{l + exp[(M - ^„a.)/T]}. By 
means of the Poisson summation formula, 13] the triple 
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sum with respect to n, a, v is transformed into another 
triple sum, = r=-oo ^vqr^ where ^^qr is expressed 
as the Fourier integral which we can calculate analytically 
by assuming /i ^ T and ujb- With these calculations 
carried out, the free energy is given as follows. 
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= *(2p^2y/^^-) cos(2p^^/c<j,) , 



(6) 
(7) 

(8) 
(9) 
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where ^'(x) = x/sinh(a;), and in eq. the coef- 
ficients (5i,Q|,Q3 take the set of indices {i,j,k) = 
(1, 2, 3), (2, 3, 1), (3, 1, 2), respectively. The summation in 
eq.(|l^) is executed further and we obtain 



QUp) = 
Qsip) - 



Stt^ p2 tan(7rpti;2/'^i) sin(7rpw3/a;i) 
Stt^ sui{Trpuj3 / UJ2) sm{TrpuJi/LU2) 



Stt^ p2 sin(7rpwi/a;3) tan(7rpa;2/'^3) 



,(11) 
. (12) 
(13) 



The free energy ilo in eq. (^ is the nonoscillating part 
with respect to the magnetic field, while (l — ili+U2+i^3 
in eq. represents the oscillating part of the total 
free energy 11. For the oscillating part, the free en- 
ergy Qi = E<tEp^i(p) vanishes because I]<t^i(p) 
Y„Qi{p) is identically zero due to eq. d). As for 
fi2, ^{p) has the oscillating factor, cos(2p7r/i/w±), in 
eq. (y) as a function of the field parameter ujb since 
uj± varies monotonically with respect to wb- Although 
the remaining factor Q^{p) in r^^(p) also oscillates by 
sin(7rpa;3/w±)sin(7rpijji/a;±) in eq. (p^, we can neglect 
this oscillation, compared with cos(2p7r/x/ti;±), since we 
can assume lo^ <^ ^ for the sufficiently large size of the 
system and also uji <^ ^ due to the assumption, ujb- 
As for ^3, fl^{p) in eq. (|^) has only a negligible oscilla- 
tion caused by sin(7rpa;i/ti;3)tan(7rpa;2/'^3) in eq. (p^), be- 
cause cos(2p7r/i/ci;3) is independent of the magnetic field. 
Consequently, fl2 is dominant in eq. (||) in terms of the 
oscillation of magnetization. Hereafter, we write the ef- 
fective part of the free energy (l as 



fin 



= E E Q2ip)'f{2p^^T/cj,) cos(2p7rM/c^,) , (14) 



where oj^ = a;2(cr). 

It should be noted, however, that eqs. (pl])-(p^ have 
spurious divergent singularities when cji/ujj is a rational 
number. To eliminate such divergences, we rewrite f2 as 
defined in eqs. (P|)-(P^, by rearranging terms, into 



n=T \ + E ^li^r , (15) 



IT — ± \p— 1 



i,q.r—l 



UJ1LU2LLI3 



(4^4)- 



.(-1)" 



4=1 



:(-i) 



(16) 



(17) 



4^-1 



-(-1)P+'' 



LOxUJ2U!'i 

r(p/i^i)-y(r/L^3) _ Y{qlu}^)-Y(Tju>3) 



[p/uJiY - {q/Ld2Y 



(18) 



where Y{x) = *(27r2Tx) cos(27r^a;), and in cq.(|7|), 
the index (j,j) takes (1,2), (2, 3) and (3,1). Equa- 
tions (p^)-([l8|) are apparently free from divergences. It 
is expected that as the temperature T rises, the to- 
tal amplitude of oscillation is rapidly suppressed by 
the exponentially decreasing factor, \E'(2p7r^r/ijjj) ~ 
{Apir'^T / uoi) exp(-2p7rr/w,). 

At sufficiently high temperatures, it is expected that 
ilo in eq. (j^) makes a dominant contribution to the total 
free energy 11, and the magnetization Mq is given by 



A^o/mb = -2 {dno/dLUB)rj.^ = -j^fi'^LOB, 



(19) 



where fiB = e/(2m) is the Bohr magneton. The mag- 
netization Mo is negative and its absolute value in- 
creases linearly with increasing lob, which is character- 
istic of the LD in a bulk system. |^ On the other hand, 
at low temperatures, by using the free energy iloff in 
eq. (p^, we obtain an oscillating part of the magneti- 
zation, M / fiB = —'2{dflcs/dujB)Tfi, which contains the 
trigonometric functions of 2ttpij,/(jJ+ and 27rp/z/w_. It 
should be noted that as the temperature T decreases, 
the amplitude of the oscillation of M becomes larger. 

In the following, we discuss how the three regions intro- 
duced by IF 1^ are modified in 3D. We regard the LD re- 
gion as the region where the amplitudes of the oscillations 
specified by p are much smaller than the confining energy 
hujs, which is unity in our notation. Hence, the criterion 
for the LD region is given by |(3^(p)|\I'(2p7r^T/w±) < 1 
for all p. Introducing the cutoff e (<C 1), we defined the 
LD region as the region where the following two inequal- 
ities hold simultaneously. 



'if{2pTT^T/uj+) < e, *(2p7r^T/tj_) < e. 



(20) 



=± p=i 



In this letter, we set e = 0.1%, but it will be necessary, 
however, to reexamine the value of e when our analytic 
classification is compared to any experiment. Noting 
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FIG. 1: Temperature-field (T-B) plane can be classified into 
three regions: the Landau diamagnetism (LD), the de Haas- 
van Alphen (dHvA) oscillation, and the mesoscopic fiuctua- 
tion (MF), in terms of the characteristic behavior of the mag- 
netization M(T,B). ujb is the cyclotron frequency and uj^ 
parametrizes the strength of the confinement. In the region 
indicated by P, the remarkable paramagnetism appears. 



^'{x)<Q and u!+>uj^, the region determined by eq. (|2C 
is equivalent to T > 7ti^+, where 7 — '^^^{e) / {2tt'^) « 0.5. 
As a result, the boundary between the LD and the oscil- 
latory region is given by T = 0.5 w+, as plotted by the 
solid line in Fig.L 

The oscillatory region T < 7W+ is divided into two 
subregions as shown in Fig.l; (A) 71^+ > T > "ftu-, 
and (B) 70;- > T. In subregion (A), the first inequal- 
ity in eq. ( pO| ) fails while the second one still holds. 
Hence, only (3^(p)vI'(27r^r/a;+)cos(27r/i/w+) contributes 
to fles- The effective free energy is given by — 
E^i Q^(p)*(2p7r2r/w+) cos{2pTTfi/uj+). On the other 
hand, both inequalities fail in subregion (B), where two 
types of oscillations regarding uj± contribute to flcs, and 
the free energy Cl^f^ is given by eq. (|l^) . As the tempera- 
ture T is lowered, the terms with larger p start contribut- 
ing to tlcff. 

For subregion (A), magnetization M derived from 
(l^g contains a trigonometric function having the pe- 
riod Lu^/fi. The dHvA region corresponds to the region 
"fuj- < T < ^uj+, characterized by the oscillation having 
the period uj-^-ffi. Particularly in a strong field region, 
since tj+ ^ ujb, magnetization M behaves as the peri- 
odic function of variable ujg^. This periodicity is typical 
of the dHvA oscillation in a bulk system. [|, | 

For subregion (B), magnetization M derived from 
contains a trigonometric function having the period w_ / fi 
in addition to In this region, lo^ and w_ are not 

so close to each other and the two modes fluctuate sep- 
arately. Consequently, the coexistence of the above two 
modes causes the rapid and irregular oscillation of mag- 
netization. This magnetic behavior corresponds to the 
MF. Hence, subregion (B) is regarded as the MF region. 



However, when the magnetic field is extremely weak 
in the MF region, the periods of uj±/fj, become much 
closer since uj± ~ 1 ± ^ujb- It is impossible to regard 
Lu± as distinguishable and therefore the separation of the 
two fluctuations breaks down in the above argument. In 
this case, the largest contribution to Cl in eq. (|l5| ) comes 
from the modes satisfying p : q : r ~ lji : 102 '■ 
namely, pjuJi - q/uj2 = g/w2 - r/tJs = r/uj^ - p/uji = 



for rip ^ r in eq. . From the behavior of characteristic 



frequency uoi 



we can see that the modes of 



p — 2j , q ^ j , r = j {j is a natural number) make a 
dominant contribution to magnetization AI . The largest 
contribution to M is evaluated as 



(21) 



Furthermore, we can estimate the upper bound U {p, q, r) 
of |f^p_g_rl^ using the triangle inequality and |i^(2^)| < 




1 for eq.(18), and it is found that the upper bound 
U{2j,j,j) decreases as U{2j,j,j) = j-'^U{2,l,l) for 
larger j. Hence, it is expected that a few modes with 
small j make a major contribution. The mode with j ~ I, 
^2.1,1, is plotted in Fig. 2(a), where we take T = 0.01 and 
N = 5000. The magnetization curve in Fig. 2(a) shows 
a remarkable paramagnetism at weak magnetic flelds be- 
low UJB ^ 0.02. The region corresponding to the param- 
agnetism is indicated by P in Fig.l. 

Next, we describe the procedure of the numerical cal- 
culation for magnetization. We consider magnetization 
as the total dipole moment of electrons with the confln- 
ing potential. From this standpoint, magnetization M 
is given by M = J^dV^ [r x J(r)], where J(r) is the 
current density distribution of the system. The calcu- 
lation of M reduces to that of J(r) which is given by 
J(r) = ReE.=(„,„,.)/(i?.)V':(r)7;f (-iV + eA)V'«(r). 
Here, f{E) denotes the Fermi distribution function, 
fiE) = {exp[{E - fi)/T] + l}-i. By using eq.(|), the 
current density distribution leads to J(r) = J{p,z)eip^ 
and 



J{p,z) 



V m ^ 



s—{n.a,iy) 



lub_p 
2 C 



(22) 



where ^ = m~2. By using eq.(E2h 
ment M is calculated analytically 
The z-component M is given by 



the total dipole mo- 
and theuM = Me^. 



M/fiB = - ^ fiEs) [a + (2n + |a| + 1) lub/coi] , (23) 

s— (n,Q,z^) 

which is also derived from M/pB = —2{dft/duJB)Tti = 
—2J2s fi^s){dEs/dujB) due to eq.(||). This coincidence 
of two derivations, of course, occurs regardless of the 
type of conflning potentials. In the numerical calcula- 
tion, we calculate the chemical potential flrst by solving 
N = '^gfiEg) for /i, and then compute magnetization 
M by using eq.(|23 
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FIG. 2: Figure 2(a) shows the contribution M2,i,i, which is 
the lowest mode {j = 1) in eq. (^l|). Figures 2(b)-2(d) show 
the magnetization curve obtained by the numerical calcula- 
tion for iV = 5000. 



Magnetization M calculated numerically is shown in 
Fig.2 for N = 5000. At low temperatures (T = 0.01, 
Figs. 2(b) and 2(c)), the magnetization oscillates over 
the entire range of magnetic fields. At extremely weak 
fields below cub 0.1 (Fig. 2(b)), the magnetization os- 
cillates smoothly with significantly large amplitude. The 
smoothness implies that not so many modes are domi- 
nant in the oscillation. Moreover, Fig. 2(b) agrees qualita- 
tively with the behavior of M2,i^i{T, ^, B) (Fig. 2(a)) with 
the chemical potential /i(T, N, B) calculated numerically 
for T = 0.01 and N = 5000. The large smooth oscilla- 
tion in Fig. 2(b) is consistent with our analysis in terms 
of free energy. At weak magnetic fields above uib 0.1, 
the large oscillation is suppressed and the magnetization 
irregularly takes a positive or negative value. This fluc- 
tuation continues until ujb 5 with the linear magne- 
tization superposed (Fig. 2(c)). At high magnetic fields 
above ~ 5, the large smooth oscillation is clear in 
Fig. 2(c). It is considered that this magnetization curve 
reflects the dHvA oscillation. On the other hand, at high 
temperatures (T = 2, Fig. 2(d)), the magnetization is lin- 
ear in lob until lob ~ 10 and then the dHvA oscillation 
increases gradually. 

Compared with the magnetization in 2D reported 
by IF, 1^ it is found that the period of oscillation is 
longer in 3D. For the 2D system, they also calculated 
the free energy analytically and its oscillating part is 
expressed as f2 = J2f=iJ2aJ2p^i (p), and fl'^{p) = 
Q1 {p)'^ {2ptt'^T / uji) sin{2pTTfi/LUi), in our notation. Here, 
we find that for 2D and 3D, these oscillating factors 
have the common phase 2p'K^,/L0i. Noting the common 
phase, when ojb is varied, the period with respect to lo^^ ^ 
A(l/ci;i), is given by A(l/a;i) = p~^fj.~^. Hence, 
serves as the scale for the period of oscillation and the 
difference between the periods can be explained by that 
between the chemical potentials. The chemical poten- 



(a) cos/co3=0.01, 77co3=0.01 (b) cog/C03=l, J'/C03=0.01 




FIG. 3: Figures 3(a)-3(d) show the current density distribu- 
tion J(p, z) for = 5000 obtained by the numerical calcula- 
tion, p is radius and z is height parallel to the magnetic field 
B. Here, ^ = l/yfmuJs and Jo = ea;3/(47rf ^). Four figures 
correspond to the characteristic regions of magnetization. 

tials /i3D,Ai2D are estimated hy N —(dflo/dfijrB and 
then /i3D ~ 6iV < fi2D ~ V2N for the same number of 
electrons. Consequently, the period of oscillation in 3D 
is longer than that in 2D. 

The current density distribution J(p, z) is computed 
by using eq.(p2[) and is shown in Fig. 3 for N — 5000. 
The positive value of the vertical axis corresponds to dia- 
magnetism. When the remarkable paramagnetism ap- 
pears (Fig. 3(a)), a paramagnetic region for J{p,z) ex- 
tends throughout the spherical quantum dot except for 
a small region near the z-axis (p = 0). In the MF region 
(Fig. 3(b)), the paramagnetic and diamagnetic regions are 
distributed irregularly. In the LD region (Fig. 3(c)), the 
current density distribution is localized near the surface 
of the quantum dot. This surface current corresponds 
to the edge current obtained by IF ^ for the 2D sys- 
tem. In the dHvA region (Fig. 3(d)), the current density 
distribution oscillates smoothly with a large amplitude. 

In conclusion, we have studied orbital magnetism in 
a 3D quantum dot with a parabolic confining potential. 
From the analysis of the free energy, we have determined 
the boundaries between the LD, dHvA, and MF regions. 
We have also computed the magnetization of the system 
and then the current density distribution, which indeed 
reflects the character of each region. The results of nu- 
merical calculation show a large paramagnetism when 
the magnetic fields are extremely weak in the MF region. 
The dHvA oscillation has a longer period in 3D than 
that in a 2D quantum dot. These results are explained 
in terms of the free energy. Furthermore, we have found 
surface current in the LD region, corresponding to edge 
current for 2D quantum dots. 
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